arXiv:1506.07500v2 [hep-ph] 3 Sep 2015 


Investigating the domain of validity of the Gubser solution to the 

Boltzmann equation 

Ulrich Heinz and Mauricio Martinez 
Department of Physics, The Ohio State University, 

Columbus, OH 43210 United States 

(Dated: September 7, 2015) 

Abstract 

We study the evolution of the one particle distribution function that solves exactly the relativistic 
Boltzmann equation within the relaxation time approximation for a conformal system undergoing 
simultaneously azimuthally symmetric transverse and boost-invariant longitudinal expansion. We 
show, for arbitrary values of the shear viscosity to entropy density ratio, that the distribution 
function can become negative in certain kinematic regions of the available phase space depending 
on the boundary conditions. For thermal equilibrium initial conditions, we determine numerically 
the physical boundary in phase space where the distribution function is always positive definite. 
The requirement of positivity of this particular exact solution restricts its domain of validity, and 
it imposes physical constraints on its applicability. 
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I. INTRODUCTION 


The dynamics and transport properties of rarefied gases and fluids are usually described 
in terms of the Boltzmann equation. The Boltzmann equation is a partial integro-differential 
equation for the distribution function f(x,p). In general this equation is solved numerically, 
and very few exact solutions are known in the literature. Nevertheless, for highly symmetric 
systems it is possible to solve this equation analytically under certain approximations for the 
collisional kernel [1], Using a relaxation time approximation (RTA) for the collisional kernel, 
the relativistic Boltzmann equation has been solved exactly for systems undergoing Bjorken 
flow [2] and Gubser flow [3, 4], The exact solution for the Bjorken flow has been useful 
to understand aspects of the isotropization and thermalization problem of a plasma formed 
by quarks and gluons (QGP) [2-10]. In contrast to Bjorken flow [11], where the system 
expands in boost invariant fashion only along one direction (the “longitudinal” direction), 
Gubser flow [12, 13] describes systems that undergo additionally simultaneous azymuthally 
symmetric expansion in the transverse directions. The solution of the Boltzmann equation 
for the Gubser flow [3, 4] was found by exploiting the *S'C>(3) ® 50(1,1) <E> Z 2 symmetry of 
the flow velocity profile [12, 13] which becomes manifest when mapping Minkowski space, 

® R, conformally onto de Sitter space times a line, dS^^R [12, 13]. In Refs. [3, 4] it was 
noticed that the resulting solutions for moments of the distribution function, such as the 
energy density or temperature, can become complex and therefore physically meaningless 
when propagating backwards in the de Sitter time (see discussion in Appendix B of Ref. [4]). 
In this work we revisit this issue and find that the unphysical behavior of the moments of 
the distribution function found in Ref. [4] is rooted in a violation of the positivity of the 
distribution function in some regions of phase space when propagating the solution of the 
Boltzmann equation for equilibrium initial conditions backward in the de Sitter time. In 
Minkowski coordinates, this translates to the distribution function becoming negative in 
certain momentum ranges at the outer edge of the spatial density profile at fixed time. 

This work is organized as follows: in Sect. II we review the procedure used in Refs. [3, 4] 
to find the exact solution of the Boltzmann equation [3, 4], In Sect. Ill we show numeri¬ 
cal results for the phase space evolution of the distribution function. Our conclusions are 
presented in Sect. IV. 
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II. THE ANALYTICAL SOLUTION OF THE RTA BOLTZMANN EQUATION 
FOR THE GUBSER EXPANSION 


In this section we briefly review the derivation of the exact solution to the RTA Boltzmann 
equation that is invariant under the group of symmetries of the Gubser flow, i.e., under the 
SO(3) q <E> SO(l, 1) (?) Z 2 group (“Gubser group”) [12, 13]. For additional technical details of 
the method discussed here we refer the reader to Ref. [4]. 

We use the following notation. The metric signature is taken to be the “mostly plus” 
convention. In Minkowski space the line element is written in Milne coordinates x M = 
(r, r, 0, q) as 

ds 2 = (y w dx ,1 dx u = —dr 2 + dr 2 + r 2 d(j) 2 + r 2 dq 2 . (1) 


where the longitudinal proper time r, the spacetime rapidity q, the transverse radius 
the azymuthal angle 0 are defined in terms of the usual cartesian coordinates (£, x, y , 

t = \Zt 2 —z 2 , q = tanh -1 

r = \Jx 2 +y 2 , 0 = arctan . 


r and 
z ) as 

( 2 ) 


The flow velocity u M is normalized as — — 1. 


A. The Gubser flow 

The dynamics of an expanding conformal fluid in Minkowski space can be understood 
in terms of a static conformal fluid defined in a particular curved space. I 11 Minkowski 
space, the Gubser flow describes a system which expands azymuthally symmetrically in the 
transverse plane and at the same time in a boost invariant manner along the longitudinal 
direction. After applying a conformal map between Minkowski space and the de Sitter space 
times a line, dS 3 (g) R, the Gubser fluid velocity becomes static in this curved space. Thus, 
it is easiest to first study the dynamics in de Sitter space, and then map the solution back 
to Minkowski space [12, 13]. 

In the de Sitter space, the Gubser symmetry can be made manifest by a suitable choice of 
coordinates. One relates the Milne coordinates = (r, r, 0, q) defined in Minkowski space 
to a coordinate system xJ l = (p, 9, 0, q) in dS;> ® R as follows [13]: One first performs a Weyl 
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rescaling of the metric: 


ds 2 —dr 2 + dr 2 + r 2 d(j) 2 


d§ =^= r z 

Next, we make the following change of variables: 

'1 - f 2 + f 2 


+ dq 2 


(3) 


p(r, r ) = — arcsinh 


9 (t, r ) = arctan 


2f 


2r 


1 + r 2 — r 2 


p G (—oo, oo) , 
0G (0,2tt), 


(4a) 

(4b) 


where r = qr and r = r/r with q being an energy scale which sets the typical transverse size 
of the system [12, 13]. Substituting Eq. (4) into Eq. (3) one finds 


ds 2 = —dp 2 + cosh 2 p (dO 2 + sin 2 6 dq i> 2 ) + dq 2 . 


(5) 


In these coordinates the metric in the curved dS 3 0 R space is then given by g /Jt , = 
diag(—1, cosh 2 p, cosh 2 p sin 2 #, 1). Eq. (5) shows that the variable p is a time-like variable. 

The Gubser flow is defined in dS 3 ® R as the unit vector = (1,0,0,0). It is the 
only time-like unit vector which is invariant under the generators of the Gubser symmetry 
group [12, 13]. The fluid velocity in Minkowski space is obtained by going back from the 
coordinates x^ in dS 3 0 R to the coordinates x M in R 3 0 R. This involves a Weyl rescaling of 
the fluid velocity components [12, 13]. One finds that in Milne coordinates the fluid velocity 
in Minkowski space is = (cosh k(t, r), sinh «(r, f), 0, 0), with the transverse flow rapidity 
[12, 13] 


k(t, r ) = tanh" 
This gives rise to the radial velocity profile 


2 rr 


1 + r 2 + r 5 


u(r,r) = tanhAc(f, f) = 


2 rr 


1 + r 2 + r 2 


( 6 ) 


(7) 


B. The exact solution of the RTA Boltzmann equation 

The invariance of a system under a particular group of transformations imposes con¬ 
straints on the number of independent variables of the distribution function f(x^,Pi). For 
the Gubser group one finds that /(£ M ;p,;) = /(p; Pq, p q ) [3, 4] where p ^ = p 2 e +p |/ sin 2 9 and 
po , p ^ and are the momenta conjugate to the coordinates 9, (j) and ? in Eq. (5). Thus, 
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the relativistic RTA Boltzmann equation in de Sitter space reduces to a one-dimensional 
ordinary differential equation in de Sitter space [3, 4], 


d_ 

dp 


f{p\ Psi,Pc) 


T{p) 

c 


f(fiPs»P<) - /eq 



( 8 ) 


where p p = p p (p,pn,p ? ) = \J'Pq/ cosh 2 p + p;? , f eq {x ) = e~ x is the Boltzmann equilibrium 
distribution function, and T = c/f re i(p) and f re i(p) = r re z(r)/r are the Weyl-rescaled (unit¬ 
less) temperatures and local relaxation time. The constant c is related to the specific shear 
viscosity of the system by c = 5p/S [5, 6, 14, 15] where r] is the shear viscosity and S is the 
entropy density. The formal solution of this equation is [1, 2, 16] 

f(p\Psi,P<) = D (p^Po)fo(po',pl,P,) + - [ dp' D{p,pl)T{pl) f eq (p p /T(p)) , (9) 

c j PO 


where D(p,p 0 ) = exp — ff^dp' T(p')/c is the damping function, and ,/o(po; Pn 5 TP) is the 
initial distribution function at po which we choose to be the Boltzmann equilibrium distri¬ 
bution, f 0 (p 0 ,Pa,P,) = fe q (p p (po)/T(p 0 )). 

The temperature T is related to the solution for / by the dynamical Landau matching 
condition, i.e., by the requirement that i(p) = £ eq {p ) ~ T A (p) where e(p) is the energy 
density computed from the non-equilibrium distribution function /(p,p^),p ? . This matching 
condition allows to rewrite Eq. (9) as [3, 4] 


T\p) = D(p,p„)U 
where 


cosh po \ 
coshp J 


f 4 (p 0 ) + - / dp'D(p yP ')n 

C Jp 0 


(S3 fv 


H(x) = ^ ( x 2 + x A 


tanh 1 (\/l — x 2 ) 


VT 


X* 


( 10 ) 


( 11 ) 


The numerical solution of (10) is uniquely determined by the initial de Sitter time p 0 , the 
initial value T 0 = T(p 0 ), and (through c) the chosen value of p/S. In order to solve for / in 
Eq. (9), we first find the temperature T(p) by iteratively solving the integral equation (10) 
as in [5, 6], and then plug this T(p) into Eq. (9) for every de Sitter time-step p, and perform 
the p' integral on the r.h.s. Since Eq. (9) is diagonal in the momentum variables p^, p f , 
the evolution of / with p can be studied separately for each point in momentum-space. In 
Ref. [3, 4] we showed that all macroscopic moments of the distribution function /, (i.e., all 
the components of the energy-momentum tensor T pu ) can be directly obtained as p'-integrals 
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over the solution T(p) from Eq. (10), with different weight functions, without even studying 
the distribution function / in Eq. (9) itself. 

In Ref. [4] (App.B) we observed for a few specific choices of p 0 that the solution of 
Eq. (10) leads to complex-valued temperatures at large enough negative values of p — p 0 . 
This observation was generic and held for any value of rj/S. The problem appeared to 
arise only for p < p Q , i.e., only in the de Sitter past; for p > p 0 the solution for T(p) was 
always real. Obviously, a complex-valued temperature is physically meaningless. In practice 
this issue can be resolved by imposing boundary conditions in the infinitely distant past 
of the de Sitter time, i.e. for p 0 —>■ —oo [4]. It was argued in Ref. [4] that the complex 
values of the temperature might be related with the violation of the positivity condition of 
the distribution function in some regions of phase space. Clearly, if this condition is not 
satisfied, the probabilistic meaning of the distribution function is lost. 

In this work we investigate the evolution of the distribution function (9) in phase space 
for arbitrary values of po, T 0 and rj/S. Results of our studies are presented in the following 
section. Before going into the discussion of our numerical results, though, we make the 
following analytical observation: Expanding the solution (9) around p 0 f° r small p — p 0 it is 
straightforward to see, for any choice of (/)^,/3 ? ), that /(p, Pq , JT) increases for p > p 0 and 
decreases for p < p 0 . The rate of increase/decrease depends on (p^,pj. From this we see 
that, as the system evolves forward in de Sitter time p, / will remain positive if it is so at 
p = po, but that / will decrease, and can for any choice of (p^,p ? ) turn negative, as we 
follow the system backward in p to sufficiently large negative values of p — p 0 . 

III. RESULTS 

In Fig. 1 and 2 we show contour plots of the distribution function /(p, p^,p f ) as a function 
of pq and p ? for fixed values of the de Sitter time p, for the specific parameter choices 
Attij/S = 3 and T 0 — 1. In Fig. 1 we impose equilibrium boundary conditions at p 0 = 0 and 
then study / at p = —4 (top panel) and p = —6 (bottom panel). In Fig. 2 we translate 
the entire problem by two units in de Sitter time p, imposing the same initial conditions at 
Po = —2 and then studying / at p = —6 (top panel) and p = —8 (bottom panel). Obviously, 
the solution is not translationally invariant in p. While in both figures we observe that the 
distribution function becomes negative at large values of p^ and small values of p ? , and that 
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FIG. 1: (Color online) Snapshots of the two-dimensional slice of the pn and evolution of the 
phase space distribution /(p,pn,p ? ) as a function of pn and p f for fixed values of p = —4 (top 
panel) and p = — 6. In both figures we consider po = 0, 47 tij/S = 3 and To = 1. 

this problem becomes more severe as \p — po\ increases, i.e. as p becomes more negative, the 
problem is clearly, for the same value of |p — po|, more serious for smaller initial values po- 
The fact that the distribution function goes negative in some region of phase space is 
independent of the value of p/S. This is seen in Fig. 3 where we show / = const, contours 
in the pn — p q plane for two different values of p/S, namely Anp/S = 3 (left panel) and 
Anp/S = 10 (right panel). In both panels we imposed equilibrium initial conditions with 
T 0 = 1 at po = 0 and plotted / at p = —6. For fixed |p — po|, as p/S increases the line 
/(p, Pq, pc;) = 0 separating physical (/ > 0) from unphysical behavior (/ < 0) is seen to move 
closer to the p q = 0 axis and towards larger pn values. In other words, for larger specific 
shear viscosity, the p-evolution that eventually drives / negative proceeds more slowly. 
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FIG. 2: (Color online) Snapshots of the two-dimensional slice of the pn and p q evolution of the 
phase space distribution f(p,pn,p^) as a function of pq and p^ for fixed values of p = — 6 (top 
panel) and p = —8. In both figures we consider po = —2, 4irp/S = 3 and To = 1. 

While the positivity condition of the phase space distribution is violated for any value 
of the temperature T 0 , the speed with which this happens as p evolves backwards from 
po depends strongly on the choice of To even if p/S is held constant. This is perhaps 
not unexpected since T also controls the scattering rate (due to the conformal relation 
Tf =const.). In Fig. 4 we present contours of the distribution function in the pQ — p q plane 
for two additional values of the initial temperature, T 0 = 2 (left panel) and T 0 = 3 (right 
panel) (in addition to the case To = 1 shown in Fig. 3), all for the same value Airp/S = 3. 
As in Fig. 3, equilibrium initial conditions are implemented at po = 0, and / is plotted for 
p = —6. As To is increased, the line / = 0 separating the physical from the unphysical region 
moves to smaller p ^ and larger p ? values, i.e. the unphysical region grows. Increasing the 




























FIG. 3: (Color online) Contour plots of f(p,pQ,p^) for fixed values po = 0 and p = — 6 but different 
values of ( 4 ^) 77 /<5 = 3 (left panel) and (4ir)ri/S = 10 (right panel). In both panels To = 1. 



Pa 



Pa 


FIG. 4: (Color online) Contour plots of f(p,pQ,p<;) for a fixed values of po = 0 and p = — 6 but 
different values of To = 2 (left panel) and To = 3 (right panel). In both panels (47 r)p/S = 3. 

initial temperature and energy density obviously speeds up the evolution towards unphysical 
behavior as de Sitter times evolves backward. 

In Fig. 5 we show this p evolution by plotting the / = 0 surface separating the physical 
from the unphysical regions in the 3-dimensional p — po, — p ? space, for initial conditions 
imposed at p 0 = 0 (left panel) and p 0 — — 1, respectively. In this figure we do not show the 
region where p ? < 0 since the distribution function is invariant under reflexions of the p q 
variable by construction [3, 4], Clearly the unphysical region only appears for p < p 0 , first 
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FIG. 5: (Color online) 3D surface defined by f{p,pn,p^) = 0 for two different initial conditions 
po = 0 (left panel) and po = — 1 (right panel). In both cases we chose (47r)p/5 = 3 and To = 1. 
The gray lines drawn over the physical boundary condition / = 0 correspond to constant values of 
p — pn — p<; over the surface. 

at large pn and later (i.e. for more negative values of p — p 0 ), also for smaller p^ values, 
and it is largest when p ? is small. As the initial time po is decreased, the unphysical region 
appears more quickly as a function of p, but covers a smaller region in pn and p q . 

The non-positive behavior of /(p,pn,p ? ) hi the de Sitter space imposes limitations of 
the validity of this particular solution when the information is mapped back to Minkowski 
space. For massless particles the distribution function is a Lorentz and Weyl scalar [4], so 
one only needs to transform the phase space coordinates (p,Po,p f ) to the corresponding 
ones in Minkowski space (p(r, t),Pq(t, r,p T ,p r ),p ? ). When transforming the momentum 
coordinates p^ in de Sitter to the corresponding ones in Minkowski space one must perform 
the associated Weyl rescaling (p M = r 2 p^) together with the covariant transformation of the 
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momentum components 


/ p p \ 

p 1 
p ( 

\l>°J 


(% % 0 0\ (r 2 v T \ 


f f 0 0 

or or 

0 0 10 

\ 0 0 0 1 / 


‘ p 
T 2 p r 
T 2 p<t> 

\rV/ 


By using the coordinate transformations (4) we obtain explicitly 


p p = T'y ( p T — v(r,r)p r ) , 


P 9 = 


T 


cosh p(r, r 


1 ip T — v(r,r)p T ) , 


y — r 2 p , 


p 2 — r 2 p' — , 


( 12 ) 


(13a) 

(13b) 

(13c) 

(13d) 


where v(r, r ) is the radial Gubser flow velocity in Minkowski space, given in Eq. (7), and 
7 = l/\/l — v 2 . Eqs. (13a) and (13b) show that up to a Jacobian factor, related to the 
Weyl-rescaling, (p p ,p e ) are related with the Minkowski-space components (p T ,p r ) by a radial 
boost with velocity v(r, r). Eqs. (13) are inverted by 


/■y 

p T = — (p p + v(t, r) cosh pp e ) , 

7 

r 


ry 

p r = — (coslipp 0 + u(r, r) , 


P0 = r 2 p* = — p* , 


Pc = P* , 


(14a) 

(14b) 

(14c) 

(14d) 


Thus, we can understand the momentum p p in the lab frame in terms of the momentum in 
de Sitter space boosted by radial flow. 

In order to develop some intuition what the break-down of the solution for the distribution 
function in de Sitter space implies in Minkowski space we make the following considerations: 

• The scale of energy q is determined by the transverse size i?j_ of the nucleus; we assume 
R± = 5 fm and thus set q = l/i?j_ = 0.04 GeV. 

• We study the system in Minkowski space at a fixed longitudinal proper time r for 
which we take f = qr = 0.3, i.e. t— 1.5 fm/c. 
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• The rotational invariance of the Gubser flow about the longitudinal axis allows us to 
choose p 0 — 0. Moreover, this lets us identify the radial component of the momentum 
in Milne coordinates with the transverse momentum of the particle. 1 


• Ps — P<; is related to the longitudinal component of the momentum in the boosted 
frame, i.e. p q = p z t — E z = w. It is an invariant under longitudinal boosts, i.e. under 
the subgroup SO(l, 1) of the Gubser group [3, 4], Thus we can evaluate p q at z = 0 
where t and r coincide and thus p q = rp z . 


For z — 0 and p^ = 0 we can write the SO(3) q invariant p ^ as 


~2 -2 , P<t> 

Pn=Pe + 


sin 

= (p e cosh 2 p(r, r)) 2 
= cosh 2 p(r, r) t 2 [y(pT ~ v(t, r) p T )] 2 , 


(15) 


where p T = a Jp\ + p\ and px = p r . 

The Minkowski-space temperature T(r, r) and its de Sitter analog T(p) are related by 
Weyl rescaling as 


T(r,r) = 


T(p(r,r)) 


(16) 


For r = 1.5 frn/c and r = 0, Eq. (4a) gives us that p ~ —1.2. When using for the 
initial condition of the temperature in de Sitter space the natural scale T(p 0 ) = 1, we 
obtain from the numerical solution of Eq. (10) that T(p = —1.2) = 0.22. Thus the 
corresponding central temperature T(r = 1.5 fm/c, r = 0) = 0.029 GeV. 


In order to visualize the distribution function in Minkowski space, we write 


f(r,r,z = 0,p T ,p z ) = f {p(r,r),p 2 n {r,r,p T ,p z ),p <; {T,p z )) . (17) 


1 In polar coordinates the momentum component p^ is 

P 4, = ~ sin (^ ~ 0p) ■ 

where pr is the transverse momentum of the particle and <j> p is the angle between the two-dimensional 
transverse momentum vector px and the horizontal axis in the transverse plane. Therefore fixing p^ = 0 
is equivalent to saying that the radial unit vector in the transverse plane is aligned with px, i.e., <j> = 4> p . 
It is straightforward to conclude that in this case p r = pr- 
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and invert the functional dependence on the right hand side. At r = 1.5 fm/c, and z = p$ = 0 
we find 


r = 5y/—0.6 sinh p — 0.91 fm , 
p z = 0.13 p q GeV, 

| 1 °- 13 7 ( J§r p + v y/ 3b +^) GeV if Pt > 7 v |Pz 
+ GeV ifp r < 7 n|p, 

where 


(18a) 

(18b) 

(18c) 


a/— 0.6 sinh p — 0.91 

0.3 — sinh p ’ , irA 

0.3 — sinh p 
cosh p 

Note that r > 0 corresponds to p < 0 at this value of r. 

In order to reconstruct the physical boundary in the r — px — p z phase space, we first 
determine numerically the set of coordinates (p,pn,ik) where / = 0 for thermal initial 
conditions at p 0 — 0, using (Att)t)/S = 3 from Fig. 5 (left panel). Once this information 
is determined in de Sitter space, it is simply mapped to Minkowski space by evaluating 
Eqs. (18) for each selected data point (p,pn,p q )- Note that the range of (p,pn,p?) covered in 
Fig. 5 limits the range in (r,pT,p z ) that we can access in this way. To extend the ( r,px,p z ) 
range requires extending the numerical results shown in Fig. 5 to a larger range in (p,Pq,P ? ) 
which is numerically costly. For this reason we show in Fig. 6 only a rather limited section 
of the surface in {r,p T ,p z ) space that separates regions of positive and negative values for 
the distribution function /. 

In Fig. 6 , the region where / < 0 lies below the surface while the physical region / > 0 
lies above it. To determine the correct value for px from (18c) we need to check the sign 
of pt — 7 ^ ; |p^|- All points shown in Fig. 6 correspond to pr > 7 vp z , i.e. to the upper sign 
in Eq. (18c). We checked numerically that the opposite sign applies in regions of large r, 
beyond the range shown in Fig. 6 . 

From Fig. 6 , we observe that in the shown range r e [2 fm, 10 fm] the distribution 
function becomes negative for sufficiently large values of pr and sufficiently small values of 
p z . As we move out to larger r values (i.e. towards the tail of the density distribution) 
the unphysical region moves towards smaller px values and thus covers a larger fraction of 
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FIG. 6: (Color online) 3D surface defined by f(r,px,p z ) = 0 for (4n)r}/S = 3 at t = 1.5 fm/c and 
q = 0.04 GeV. The gray lines drawn over the physical boundary condition / = 0 correspond to 
constant values of r — px — Pz over the surface. See text for further details regarding the initial 
conditions. 

momentum space. For larger r, the solution with thermal boundary conditions at p 0 = 0 
thus remains physical only for very small px values combined with sufficiently large \p z \. At 
fixed r, the critical surface separating / > 0 from / < 0 is almost linear in px and p z and 
approximately given by \p c J lt \ = a{r){px — Pt q (^)) where px 0 ( r ) is the li nc where the critical 
surface intersects the p z = 0 plane. We fold a(2 fm) « 1.47 and a(10 fm) « 2.08. These 
values depend on the initial choice of r and all the other chosen parameters (T 0 = T(p 0 ), 
p/S and q ). For p z > |pf nt | the solution is physical (i.e. / > 0) while for \p z \ < |p^ rit | it is 
unphysical (/ < 0). 

IV. CONCLUSIONS 

In this work we have studied the dynamics in phase space of the recently found exact 
solution to the conformal RTA Boltzmann equation which undergoes Gubser expansion [3, 4]. 
We determined the distribution function by first solving equation (10) for the temperature 
T(p) and then evaluating Eq. (9) point by point as a function of the variables (p, PmP?) that 
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define the phase space of the system in de Sitter space. We assumed thermal equilibrium 
boundary conditions at de Sitter time p 0 . 

We observe that this exact solution becomes negative, and therefore physically meaning¬ 
less, for sufficiently large negative values for p — p 0 . This proves the conjecture in [4] that 
the unphysical behaviour of the moments of the distribution function is likely related with 
the violation of positivity of the distribution function in certain phase space regions. 

The non-physical behaviour of the distribution function depends strongly on the initial 
value p 0 . If the evolution of /(p,pn,p ? ) starts at p 0 = —oo the system is always evolving in 
the forward direction in the de Sitter time p and the distribution function remains always 
positive. When imposing thermal equilibrium initial conditions at some finite value po, the 
system can evolve not only forward but also backward in de Sitter time p. In this case the 
distribution function /(p,pn,p ? ) is always positive definite for p — p 0 > 0 but it loses its 
probabilistic and physical meaning in some regions of the phase space when p — p 0 < 0. 
Generically, the regions where the distribution function is negative occur at small p f and 
large p^, qualitatively independent of the value of the shear viscosity over entropy density 
ratio r)/S and the initial temperature To. The generic shape of the surface separating the 
physical (/ > 0) from the unphysical region is shown in de Sitter coordinates in Fig. 5 and 
in Minkowski coordinates in Fig. 6. In Minkowski space, the exact solution for / becomes 
unphysical at large r, large Pt and small p z \ this corresponds to large negative values of 
p — p 0 , large p^ and small p f . The example studied in Sect. Ill shows that when choosing 
values for the parameters q, T 0 and rj/S that are natural for heavy ion collisions, problems 
of non-positivity of / arise already at moderately small values of r and px when p z = 0. 
This renders the analytical solution found in [3, 4], which assumes local thermal equilibrium 
on a surface po, unsuitable for heavy-ion phenomenology. 

This does not mean, however, that this exact solution cannot be used to test different 
hydrodynamic approximation methods, as done in Refs. [3, 4], There is no problem with 
such tests as long as the comparison is performed (either in de Sitter or Minkowski space) 
in the physically allowed region where / > 0. An alternative approach which guarantees 
always the positivity of the exact solution is to fix the initial condition at sufficiently large 
negative p 0 values such that the region p > p 0 includes all of the interesting range in r and 
r covered by the evolution of a heavy-ion collision. In this case, with a suitable choice for 
the initial form /o(po>Pn,p f ) that is presumably not of equilibrium form, the exact solution 
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of Eq. (9) might become relevant for heavy-ion phenomenology. 
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